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solution  of  conjugate  heat  transfer  including  convection  in  a fluid 
medium  and  conduction  in  the  phase  change  material,  a problem 
applicable  to  tho  analysis  of  latent  heat  thermal  energy  storage 


INTRODUCTION 


smoothly  dispersed.  Probably  most  intriguing  of  all  is  the  phase 
change  in  alloys.  In  the  alloys,  the  phase-change  temperature  is  not 
distinct.  There  is  a crystal  line- like  structure  at  the  interface  which 


conveniently  with  computational  approximations.  This  is  not  so. 


Cinpeyron  (see  Rubinstein,  1971)  first  st' 


disciplines.  It  is  t 


2.1.1  Exact  solutions 


only  a few  terns  are  thus  sufficient  for  an  accurate  solution. 


differential  equations.  Lightfoot  (1929)  developed  e 


function.  Since  Lightfoot  ui 


obtained  in  bis  work  can  be  solved  by  a similarity  transformation, 
functions  (Crank,  1984).  Thus,  Lightfoot's  intcgrodiffcrential 

leads  to  the  use  of  a moving  heat  sink  for  melting.  Lightfoot's  method 
can  thus  be  taken  as  a source-and-sink  method.  Chuang  and  Ssckcly 


1992).  In  these  efforts. 


(llsieh  and  Choi,  1992a).  In  th 


n important  feature  o 


Furthermore,  in  the  source-nnd-sink  method,  the  interface  flux 


(1974)  in  the  solution 


a fixed  domain.  However,  since  the  conditions  imposed  on  the 


introduced  a fictitious  flux  condition  at  the  fixed  boundary.  This 


Stefan  problems  can  also  be  solved  by  a heat  balance  integral 


commonly  known  as  the  method  of  weighted  residuals,  the  heat  integral 


method,  the  governing  equation  is  changed  to  an 


appropriateness  of  this  function  in  the  solution.  Goodman  and  Shea 

change  probcms  with  two  phase  regions. 

It  has  been  widely  accepted  that,  in  the  heat  balance  integral 
method,  using  a higher-order  polynomial  for  temperature  in  the  trial 

solution  (Goodman,  1964;  tangford,  1973).  However,  a spatial  sub- 


analysis of  solidification  around  a cylindrical  pipe  by  a Runge-Kutta 


(19G9)  have  demonstrated  that  the  quasi -stationary  approximation  is 
actually  the  zeroth  order  solution  of  a regular  perturbation,  which 


and  Jiji  (1977)  have  shown  thi 


prohibitively  tedious,  the  method  is 


reduced  to  an  initial-value  problem.  In  a similar  fashion,  a Laplace 


analytical  results  were  in  reasonable  agreement  with  the  numerical 


Stefan  problems  con  also  be  solved  by  using  an  embedding  method  os 
shown  by  Sikarski  and  Bolcy  (1965).  Following  a totally  different 
approach.  Rathjcn  and  Jiji  (1971)  were  able  to  combine  Lightfoot's 

condition  (Patel.  196S).  They  developed  an  analytical  solution  for 


Solution,  of  the  Stcf.n  Problem 


In  the  solution  of  the  phase  change  by  a classical  formulation, 

interface  position  any  not  fall  on  the  grid  lines  initially  set  up  in 
the  solution,  iteration  is  usually  necessary  which  lends  to  complexity 


transformation  of  variables.  In  these  methods,  the  interface  position 


is  always  tied  to  a grid  line  or  is  immobilized  in  the  transformed 


Wilson  el  al.  (1978) 


popular  methods  have  been  provided  in  works  by 
Crank  (1981). 


the  previous  time  step  are  mapped  onto  the  new  grid  by  interpolation. 
These  methods  are  often  referred  to  as  front-tracking  methods  and  have 
been  used  extensively  in  the  numerical  solution  of  phase-change 


solution  of  solidification  in  continuous  casting.  Additional  numerical 
Saitoh  (1978);  Pham  (1989)  and  Liang  (1993),  among  others 


2.2.2  Weak  formulation 

analysis,  the  weak  formulation  can  be  usee 


(Elliot 


survey  of  the  recent  developaeuts 
given  in  Turnon  end  Namburu  (1980) 


arbitrary  soootliing  function  of  the  enthalpy  (Goodrich,  1978)  and 
suppression  of  the  oscillation  in  the  solution  (Vollor  and  Cross, 


1983;  Phao,  1985,  1986).  Vhile  these  difficulties  can  be  resolved  by 


Phase-change  problems  have  recently  been  solved  by  a boundary 


T(?.‘,)=T0ir) 


T'(S,t)=T.(S,t)  = T„ 
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Eq.(3.10), 


for  [BC]  in  the  second 


Boundary  condition.* 

Expressions  for  [BC]  in  Eq.(3.10) 

T(  i,,l)  = rt(l( 

T^(r)flC(t^|x„r) 

UTM  fe(<) 

"‘(rlC(.,f|i..i) 
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+ £f  ^G(^'\SM,r)dr  (3.14) 
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superposition  method  with  o lieu  time  mnrehing  scheme  is  developed. 
Before  addressing  this  method,  the  validity  of  the  superposition 


established. 


t follows,  it 


and  equal  property  values  in  different  phase  regions.  Thus,  the 
governing  equations  are  all  linear.  The  initial  and  boundary 
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3.2.2  Solution  method 


equation  is  given 


* $=£  &*©•  ien  C3-28) 

where  zd  is  the  Sturm-Mouvi 1 le  weighting  function,  given  as  follows: 


The  initial,  boundary,  and  interface  teaperature  and  flux  conditions 


**f-**M0- 


+ (3.36) 

In  this  problem,  Eq.(3.26)  is  satisfied  implicitly  because  of  the 


Correspondingly,  the  V(s,l)  problem  satisfies  the  following 


Sr)-  *="••  l€tI  (3M> 


(3.38) 


analytical  solutions 


Following  Hsieh  and  Choi  (1092a),  the  interface  condition  (3.40)  is 
incorporated  into  the  governing  equation  (3.36)  to  form  an  equivalent 


ix.,  sen  (3.4i) 


(sec  Eqs.(3.37)  and  (3.3S)),  the  V problem  can  be  solved  as 

V(M)=*j  T^V(r)C(s,(|S(r),r)dr  (3.42) 

Application  of  the  interface  condition,  Eq.(3.39),  gives  a relation 


rm-„(S,„=y^ 


S,'(r)C(S(l),l|.S(r),r)rfr  (3.43) 


located  at  a grid  point,  U(S,I)  can 


time  step  after  phase  change,  the  T-problcm  is  decomposed  i 
problem  and  a V-problcm.  The  {/-problem  satisfies  t 
initial  and  boundary  conditions  (see  Eqs.(3.33)  through  (3.35)),  vhile 

together  with  interface  conditions  (soc  Eqs.(3.36)  through  (3. 40)).  In 


(3.36)  and  (3.40),  a reduction 


(1)  Start  from  an  initial  condition  T0(x)  and  set  n = l. 


(4)  Calculate  l^x,!"*1)  from  Eq.(3.42), 


(5)  Find  T(.,l"*>)  ={/(., 1"*')  + F(.,I"«), 

(6)  Sot  T0(«)  = T(x,l"*1),  n = n + l,  and  repeat  (2)  to  (6). 


integral  in  V(x,l) 


|.(3.42))  i< 


Solutions  of  One  dimensional  Stefan  Probli 


The  analysis  given  in  the  preceding  section  will  now  be  used  to 
solve  one  dimensional  Stefan-problem  examples  in  Cartesian,  cylin- 
drical, and  spherical  systems.  For  the  examples  in  Cartesian  system, 
both  semi-infinite  and  finite  domain  (plan  slab)  will  be  solved, 
whereas  in  the  cylindrical  and  spherical  systems,  only  finite  domain 

the  governing  equation,  initial,  boundary,  and  interface  conditions 
have  been  given  as  Eqs.(3.28)  through  (3.32).  The  V and  V problems  are 
formulated  as  Eqs.(3.33)  through  (3.40).  The  V problem  has  been  solved 
using  the  equivalent  problem  (3.41)  with  the  result  given  as  (3.42). 
The  interface  position  is  to  be  solved  using  (3.43).  Solutions  for 


e governing  equation  for  the  (/-problom  in  a Cartesian  system 
ibtained  from  Eq.(3.33)  by  setting  rf  equal  to  zero.  To  ensure 
1 stability  for  all  time  steps,  an  implicit  finite  difference 


-tjfpl  + ( -O.f'+l  = V J (3.43a) 

-frj'jtl  = &u"-,  + ( 1 -»»)!/;  + (3.45b) 


superscript  n refers 
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The  governing  equation  for  the  (/-problem  in  a cylindrical  system 
can  be  obtained  from  Eq.(3.33)  by  setting  d equal  to  unity.  An  implicit 
finite  difference  solution  for  this  problem  can  be  derived  as 


+ ( 1*28,)('J*‘  - (»,*(>))(/;;;  = (/;,  ifr^O  (3.S6) 


Here,  = oAI/(2zAr)  i 


•’ilr 


«'(*.'.)  =*j  ^3S»(r)C(«,l1|S(r,,,)*  (3.63) 

(3.64) 


2(j;+k2) 


0ncot(0nR9)  + A"  = I 


(3.67) 


)><ll  S/r),r)ir 

Tm-V{S„l,)  = ^^(r)C(SJ,l1|S,(rl,r)ir 

(3.72) 


— ^-isJ(r)C(S„.ll  |S,(r),r)dr 


.5  Stability  analysis 


stability 


r"*i  = (/n*l+i"1*i 


(3.74) 


(3.76) 


the  solution  of  U is  obtained  by  using  a finite  difference  method. 


if  the  simple  implicit  method  (3.46a) 


(3.45b)  is  used  (Anderson 


The  disturbance  <2  must  satisfy  the  following  equatic 


i Ti  = i?  * i*3i 


■»=*j  Sc(».«,|iM.r)* 


As  discussed  previously,  for  small  Al,  rfq/dr  can  be  treated  as 


v'(S,i,)=rm-i;(s,i1) 

‘’U= 


(3.84) 


(3.85) 
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Position,  x/H 


Figure 


«-0*0 


TfaS,0=T,(x,S,t)  = Tm 


Figure  3.2  System  analyzed  in  a two-dimensional 


r(i,i/,i)= <'(*,»,!)  +v(x,v.*> 
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I'jlx.s.o  = v,(x,s.o  = r ,„-</<> 


will  be  necessary  if  this  position 
depends  greatly  on  the  fitness  of  1 


directly  for  S.  The  difficulty  lies  in 


differential  i 


f V with  * in  (3.100)  is  changed  to  a finite  difference 
fj‘V(z,ll,t)  y(i*Ar,v,l)  - 2y(r,y,l)  -f 


Then,  by  indexing  i position  with  subscript  i so  that  l'(r,y,l)  is  changed 
to  Vf(p,<)  and  K(z± Ax.ir,!)  to  Kl±1(»,l).  Eq.(3.100)  can  be  re-written  as 


iflVjy.i)  a‘v,(u,i)  vmM-wM+v,-r 


-h  % 


"H*" 
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<o<'<‘l  (3.111) 


(3.112) 


Eq. (3.100)  con  be  used  to  re-urite  (3.113)  os 

> ' 23f4<‘,-S.)  +' K*(»-S.-l)^i  (3.114) 


',i(#.<l)=IJ  ^ 01(»,l,|S„r)dr 


('-<„)-p!C,(ff,<,|S.(l„),r)dr.  « = i.,± 


a first-order  approximation  is  sufficiently 


*,(»'.  t)  = ».(y',  I 


(3.117) 
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Similarly,  the  second-order  approximation  of  results  in  the 


■>-*] 


•jJ  Cj(t,l1|SJ,r)dr+ll'J 


(3.13d) 


“’«=j  ('-'o)^ra!Oj(t.<ilS«|l„),r)dr. 


= s7  + ^(r3»)+!fe  % <(r-S(t,<|)  (3.137) 


’ ' A? 


■)=ij  7F*.W 


>'iCr,l1)  =|j  ^S,(r)C,(r,lI|S„r)dr 


;S.('<,)C1(r,l,|S.|l0|,r)dr 


Solutions  in  three  dimensional  Cartesian  systems 


S = S(*.  »,0 


[1+(s)J  + (ff),][^-w]  “HT"  3?  (3.144) 


T(’.U.t)  = U(i,»,.-,l)  + V'(i, (3.145) 


+ -2Pi,i,*+Pi,i,*-I)  (3. 146) 


i si'  = -^r+b  gf  0 


VU(..«.{|^ 


*a£? 


U(S,I)  which 


in  the  integrand  will  not  be  magnified  but  will  be  suppressed.  On  the 

the  ^-problem  through  the  superposition  principle.  However,  this 
disturbance  can  still  be  suppressed  by  the  ADI  method  in  the  solution 
for  the  temperature  in  the  succeeding  time  step.  It  is  thus  expected 


dimensional  problem  (See  Chapter  5)  will  be  provided  later  for  tests. 


applied  directly  to  the  phase  change  material.  In 
change  material  may  be  heated  or  cooled  by  a flu: 
itself,  a conjugate  heat  transfer  problem  tht 


Two  conjugate  heat  transfer  problems  will  be  solved  in  Sections 


a straight  parallel  c 
first.  Results  for 


Problem  formulate 


A thermal  energy  storage  (TBS)  unit  is  sketched  in  Figure  3.5. 
This  unit  can  be  modeled  for  analysis  as  shown  in  Figure  3.6.  It 
consists  of  a straight  parallel  channel  surrounded  by  PCM.  Heat 
transfer  fluid  (HTF)  is  forced  to  flow  through  the  channel  and  it  ex- 
changes heat  with  the  PCH.  It  is  a time-dependent  phase-change  problem 


accomplished  by  solving  the  Navier-Stokes  equations  and  the  energy 
equation.  However,  a direct  numerical  simulation  of  the  Navier-Stokes 
equations  is  time  consuming.  In  addition,  for  a conjugate  problem,  an 
iteration  is  necessary  to  satisfy  the  boundary  condition  at  the 
fluid/PCH  interface,  which  further  adds  to  the  complexity  of  the 
solution.  An  alternative  approach  is  taken  that  uses  empirical 
correlations  for  the  treatment  of  the  heat  transfer  at  the  fluid  side. 
This  permits  the  use  of  convective  heat  transfer  coefficient.  While 
the  accuracy  of  the  method  depends  on  the  accuracy  of  the  empirical 


mj  i ,■”> 1 hiTt-y 
7r*v^r+  pi',j 


t transfer  coefficient  A is  taken  from  the  literature.  I 


' condition  (Kays  and  Crawford,  1980).  There  is  a h 
il  equation  developed  by  Gniel inski  (1976)  as 


= [0.79/n(  WeD)  — 1.64  J-2 


e hydraulic  diameter  of  the  channel.  It  is  noted  that 


straight  channel  by  redefining  t 
Dh  as  shown  in  (3.160). 


2 Solution  method 


problem,  the  heat  flux  at  the  channel  wall  and 

by  Barozzi  and  Pagliarini  (1985),  even  for  a 
conjugate  problem  without  phase-change. 


instabilities  may 

(3.154)  is  re-written  in  ter 
grid  point  above  the  boundary 


yp  (Figure  3.6).  Thus, 


*?'/ 1 ,,eTl , TJ-T i _n 


E, = |^,[T((0.r)-T,(//,r)]dr 


Hti* 


+J[-i,s(*,oy* 

b*l«»cc  ip  (3.187)  is  fa 


Figure  3.8  Control  volume  for  fluid  analysis 


Figure  3.9  Boundary  conditions  for  PCM/I1TF  analysis 


k-i 


(3.174) 


Mr.-r/,= 

In  thin  effort,  the  convective  heat  transfer  coefficient  A for  the 
packed  aphcres  is  obtained  by  curve-fitting  the  experimental  data  of 


t transfer  coefficient  c 


Solution  of  the  conjugate  problem  for  fluid  flow  in  a packed 
sphere  bed  follous  the  same  line  described  earlier  for  channel  flow. 


is  negligibly  small,  so  that  the  heat  transfer  from  the  fluid  to  the 
sphere  can  be  expressed  in  terms  of  the  temperature  difference  between 


For  the  spherical  PCM,  the 


’-i  + ^T^-jy+T^-Ty  (3.179) 


(sec  Figure  3.9).  Equatioi 


T,  in  Eqs.  (3. 178). (3. 180)  is  replaced  by  Tm  and  r,  is  changed  to  5. 


■* 


(3.1 


represents  the  tine  rate  of  change  of  the  internal  energy  inside  the 
control  voluao.  For  the  problem  under  investigation,  there  are  energy 
changes  in  the  PCK,  the  shell  of  the  PCM  sphere,  and  the  fluid. 


«f£T/<0,r)-7yff,r»»r 


+ psK?|  Prtfr^O-T. Jfa 


f --'"-Koir  of  ^ph^e 


e cncrg>'  associated  w 


the  tank.  It  is  expected 
the  problem,  the  left  ha 


accurate  numerical  solution  of 


(3.185)  should  be  equal  to  the 


used  to  check  the  accuracy  of  a numerical  experiment  to  be  presented 


CHAPTER  4 

TVO  DIMENSIONAL  PHASE-CHANCE  EXPERIMENT 


In  this  chaptor,  an  experiment  is  described  for  the  purpose  of 
validating  the  numerical  solutions  described  in  the  previous  chapter. 
The  experimental  setup,  method  of  measurements,  and  data  reduction  are 


is  chapter.  Since  there 


r validation  of  the  numerical  solution. 


it  has  been  veil  established  in  phase-change  studies  that  the 
interface  position  provides  the  best  indication  for  the  solution 

if  the  interface  position  is  conspicuously  in  error.  This  is  due  to 
the  fact  that,  for  a phase-change  problem,  the  boundary  and  interface 


facilitate  visualization  of 


coll  of  the  configuration  of  a rectangular  parallelepiped  is 
constructed.  Paraffin  wax  is  chosen  as  the  phase  change  material 

point  (Lane,  1988).  It  also  possesses  the  desired  characteristic  of 
becoming  translucent  to  visible  light  upon  phase  change  from  solid  to 
liquid  state.  The  interface  position  can  thus  be  measured  by 

the  interface  position  on  photographic  films  by  means  of  a camera.  The 
thermophysical  properties  of  the  paraffin  wax  are  summarized  in  Table 
4.1  (sonree:  Lane,  1988),  in  which,  the  melting  point  listed  is  the 
result  of  experimental  measurements  that  are  specifically  performed 


Conductivity 


(«■) 

(w/tf) 

(kJ/t,  K ) 
(W/m  K) 
W'"3) 


is  designed  so  that  the  boundary  condition  can  be  controlled 
precisely.  Provisions  are  thus  made  to  simulate  the  prescribed 
temperature  at  the  heated  side  of  the  wax  while  keeping  the  other 


expansion  of  the  wax  upon  phase  change.  Another  provision  is  thus  made 


c gap,  providing  justification  of  use  of  a two  dimensional  t 
»nsfcr  analysis  for  the  cell  (more  justification  to  follow).  To  1 
j phase  change  material,  a heater  (/)  is  provided  at  the  top  of  ti 


Figure  4.2  Assembly  drawing  of 


specially 


heating  elements  (6)  which  are  individually  wii 


le  direction  while  maintaining  a gradual  temperature  variation 
other  direction.  Coppcr-constantan  thermocouples  (d)  arc 


holes  provided  at  the  flnngcs  of 


holes  to  the  small  recess  formed  between  the  flange  and  the  wooden 


Paraffin  wax  has  been  used  ns  the  phase  change 
poured  into  the  cell  through  the  hole  at  the  base  of 
Figure  4.3).  The  cell  is  placed  in  an  upright  position 


facing  the  light  source  (/  in  Figure 


Figure  4.4  Exploded 


providing  grid  lines  to 


Figure  4.1)  is  used  to  record  the  images  of  the  vox  during  phase 


interface  positf 
the  time  delay  t 


unavoidable  in  recording  temperatures  from  one 


4.3  Sample  Results  and  Experimental  Uncertainty 
Sample  pictures  taken  for  the  interface  positions  are  shovo  in 
Figure  4.5.  The  interface  positions  are  clearly  visible.  Also,  because 
of  the  temperature  imposed  on  the  wax  being  nonuniform,  there  are 

the  grid  paper  used  has  a scale  of  five  divisions  per  centimeter. 
Repeated  tests  of  reading  the  interface  position  over  the  background 
grids  provides  an  uncertainty  of  position  reading  to  be  about  0.5mm. 


CHAPTER  5 

RESULTS  AND  DISCUSSION 

In  this  chapter,  numerical  experiments  will  be  provided  for  the 
validation  of  the  solution  methods  described  in  Chapter  3.  Numerical 
solutions  of  one  dimensional  Stefan  problems  will  be  discussed  first. 

the  results  published  in  the  literature.  Solutions  of  two  dimensional 


Dimensional  Problems 


compared  with  analytical  s 
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readily  available.  Also,  the  solution 
,ite  domain,  Eqs.(3.62)-(3.55),  if  the 
for  solution.  As  shown  in  Figure  5.1, 


Figure  5.2  shows  the  interface  positions  versus  time  for 


s.  corresponding  to  the  Fourier-modulns  (oAl/Ar1) 


Figure  5.1  Temperature  profiles 


superposition  method 
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Case  2 (Table  5.1)  provides  a more  thorough  ; 
superposition  method  based  on  the  development  of  tl 

and  flux  conditions  are  imposed  on  a semi-infinite  medium  (aluminum), 
which  is  initially  subcooled  (see  Table  5.1).  A two-phase  melting  thus 
occurs  and  the  numerical  results  are  compared  with  the  SSH  results 
reported  by  llsich  and  Choi  (1992a)  and  Choi  (1992). 
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perspective,  the  ■odium  io  a finite  domain  with  back  side  (*=ff) 
insulated  nay  be  considered  to  be  pre-heated  prior  to  melting.  Thus, 
for  the  final  stage  of  melting,  the  medium  is  fully  pre-heated  to  its 
melting  point;  the  problem  reduces  to  a one-phase  melting  in  a Stefan- 


s represents  1 


!>  a sinusoidal 


(x-fl)  are  identical.  Clearly,  because  of  the  boundary  conditions 
imposed,  the  interface  position  rapidly  reaches  a quasi-steady  state. 
The  solid  curve  fluctuates  about  the  dashed  curve,  vhich  eventually 


there  is  approximately  a 10  degree  angle  of  phuse  delay.  The  interface 
oscillates  at  the  same  frequency  as  the  imposed  temperature  cycle.  The 
long-time  cyclic  solution  presented  here  is  believed  to  be  the  first 
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5.1.5  Melting  of  a sphere 

In  all  the  tests  performed  above,  the  phase  change  has  been 
considered  in  a Cartesian  systca.  Case  8 deals  with  a one  dimensional 
phase  change  in  a sphere  which  is  initially  at  its  melting  point  but 
in  solid  state.  The  surface  is  suddenly  heated  to  a constant 
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in  Table  5.1.  There  is  no  exact  solution  for  this  problem.  Hill  (1987) 
solved  this  problem  by  an  integral  method,  and  derived  the  lover  and 
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following  conclusions  « 


i -dimensional  phase-change  problems. 


technique.  The  validity  of  the  superposition  method  for 
the  solution  of  Stefan  problems  with  constant  and  equal  properties  in 
different  phases  has  been  rigorously  analysed.  In  the  superposition 
method,  the  solution  of  the  Stefan  problem  was  superimposed  by  a 
finite  difference  method,  which  was  used  to  calculate  the  temperature 
and  boundary  conditions,  and  an 
d to  track  the  interface  position 
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and  to  update  the  temperature  due 

combines  an  analytical  solution 
difference  method  in  anoi 
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Hating  interfaces  resulting  from  time-variant 
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in  the  diecharging/recharging  processes  of  a latent  heat  TES  unit  wore 
successfully  calculated.  Coalescing  front,  in  . tvo  dimensional 
melting  and  freezing  problem  with  multiple  phase  regions  were  also 
investigated. 

5.  A tvo  dimensional  phase-change  experiment  vas  conducted  for 
testing  the  accuracy  of  the  short-time  solutions  of  the  numerical 

long-time  solutions.  In  both  cases,  the  interface  positions  were  used 
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experimental  mi 
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e extended  to  the  solution  of  three  dimensional  phase-change  problems 
uch  as  laser  cutting,  drilling,  and  surface  hardening. 

2.  The  properties  of  the  PCM  in  the  present  analysis  have  been 
onsidcred  as  constant  and  equal  in  different  phase  regions.  In  the 
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8 completes  the  derivation  of  Eq.(3.52). 

For  a fin*  condition  at  * = 0,  the  Croon’s  function  is 


[see  Eq.(3.51)) 


G(*,<  l*Vr)  = C„+  53^t'“3"|,'')co«(a„i)co«(^„i' 


Substituting  (A. 9)  into  (A.l)  yields  the  following  for  solution 

= Cjffi&l  + |t£c„((l)co.(tfnx)  (A.  10) 

C„(l,)=|  ^{r°ti('‘-,,co.(t„S,r))}dT  (A. 11) 

By  following  the  same  line  as  derived  in  (A.5)-(A.6),  Eq.(A.ll)  can  be 
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derivation  of  Eq.(3.53). 


APPENDIX  B 

DERIVATION  OF  EQUATION  (3.68) 


|S(r),r)dr 


(B.l) 


surface  is  given  by  (see  Eq.(3.64)) 


Substituting  (B.2)  into  (B.l)  yields 


(B-3) 


(B.4) 


S(r)=S(l„)  + „( r — 1„)  , 


(B.S) 


| «"»in(ir)dr  = . sin(tz)  - 6 . „.(»*)] 

- “ ‘V-fte)  - 2ot  • coa(fa)j  (B.7) 

cn,  Eq.(B.S)  can  bo  integrated  and  rearranged  as  follovs 

D"<'1  ’ = (“S)‘  + 

" [W*)5' 'I  (dn»)1|,“S")!',a"'11]  [*'"Wn«U)  * 

-2o^.[co,((InS(l|))-,-fl5a,„.(SnS(l0|)|  (B.8) 

s completes  the  derivation  of  Eq.(3.68). 


FORTRAN  PROGRAM 


PROBLEMS 


APPENDIX  C 

FOR  TWO  DIMENSIONAL  PHASE  CHANGE 


REAO( IREAD. 1000)  TITLE 
READ(IREAD,*) 

read(iread’.) 

READ) IREAD, 1000)  FREST 
READ( IREAD, •) 

readJiread!.)  DTT,TOUT,TSTOP,KREST 
NXY  =ZRAX(NX,NV) 


CALL  I N I T( NX , NY , SBC, I BC , FREST , KREST , TSTART 

CALL  PR0CES( NX , NY , NXY , KBC , I BC , ITFS , TSTART , 
DUM2(  I ) , DUM2(NXY+1 ) ) 

0PEN( I OUT,  FILE=FOUT) 

CALL  OUTPUT(NX,NY,TSTOP,Y,T) 

CALL  DUMP(NX,NY,FMIMP.TSTOP,SO,S,Y,T,TV) 
FORMAT(AAO) 


IF(NX.EQ.l)  THEN 
DX  = XDIH/FL0AT(NX-1) 

DY  = YDIM/FL0AT(NY-1) 

ALPHA  = CK/(IUIO.CP) 

1F(KREST.EQ. 1 ) THEN 

OPEN(  I REST , FI  LEeFREST , FORM; ' UNFORMATTED  ’ , STATtlS= 
CALL  RESTART(NX,NY,TS,SO,S,Y.T,TW) 

CLOSE(IREST) 

PRINT., ’===  RESTART  OK  ! ======’ 


PUT  BOUNDARY  CONDITIONS. 


NBC  = lABS(KBC) 

IF(NBC.GT.HAXBC)  STOP  ’KBC  Is  too  large’ 
READ(2|.)  (XBC(I ) ,1=1 ,NBC) 

READ(2!.)  (TBC(I ,2),I=1 ,NBC) 


CONTINUE 

CONTINUE 

IF(TS.LE.TIB2)  COTO  500 

READ(2..,END=500)  TIB2 
READ(2,.)  (TBC(I ,2),I=1 ,NBC) 


END1 IF=  TBC(I,1)  * (TBC( I ,2)-TBC( 1,1) ).(TIME-TI1 
CONTINUE 

Y2(l)  =’-0.5 


p'G  = SICCY2(I> 
u(i)  ^ i 

CONTINUE 


-1))/(XBC(I.1)-XBC(I-1) 


6..((T8(U1)-TB(I))/(XBC(H1)-XBC(I)) 
(TB( 1 )-TB( 1-1 ))/(XBC( I)-XBC( 1-1 )))/ 
(XBC(Itl)-XBC(I-I))-SIG.U(l-l))/P 


UN  — (3./(XBC(NBC)-XBC(NBC-I))).(YPI 
1 (XBC(NBC)-XBC(NBC-I ))) 

Y2(NBC)  = (UN-QN.U(NBC-1))/(QN.Y2(NI 


l-(TB(NBC)-TB(NBC-I))/ 


IF(XX.GE.X0C(J-1).AND.XX.LE.XBC(J))  THEN 
II  = XBCfJ)-XBC(J-l) 

A = (XBC(J)-XX)/H 
B = (XX-XBC(J-1))/H 
TERM I = A.TB(J- I ).B.TB(J) 

TERM2  = ((A..3-A).Y2(J-1)+(B..3-B).Y2(J)).(II.II 


TV(I)  = TO 

CONTINUE 

RETURN 


SUBROUTINE  RESTART(NX,NY,TS,SO,S,y‘t,TW) 

CO»RONI/IO/I^AriOOTNlS^D^p'2)'S°(NX'2,'TV<NX) 

READ( IREST)  ((Y(I , J) ,T( I , J) , 1=1 ,NX) , J=1 ,NY) 

READ(IREST)  (S(l,l ) ,S(I ,2),S0(I , 1 ) ,SO(I ,2),TV(I) ,I=1,NX) 
RETURN 


E DUMP(NX,NY,FDUMP.TS,SO,S,Y,T,TV) 

0/IKAriOWNIRS’lDUNX'2)’S°<NX'2)'TV 


RM= ’ UNFORMATTED ' ) 


« 


C========  SUBROUTINE  PI 


SUBROUTINE  P 


OCES(NX,NY,N 
IT.IICOV.TINF 

COMMON  /BC/XBC(MAXBC),TBC 

COMMON  /ALL/XDIM, YDIM.DX, DY , 
COMMON  /SS/TSI1 ,TS12,TS2I ,TS! 


TIME  = TSTART 


TV1(I)  = TV(I) 
ITFS(I,1)  = 0 
ITFS(I ,2)  = 0 


°UIF(Drr.LT.O.  AND.  NCY.LE.10)  THEN 
DT  = ABS(DTT).(0.1.NCV)..Q 

^DT^=  ABS(DTT) 


WUTB(.,1000)  S(I,l),S(I,2),T(I,l),(i-|).dx 
CONTINUE  ' ' 1 


IP(IP.Eq.O)  VRITE(30,1000)  ( 


■DX,S(I,1),S(I, 


1) .S(NX.l) 

2)  ,S(NX,2) 


CONTINUE 

IF(TIME.LE.TIB2)  GOTO  500 

TBC(I,1)  = TBC(I,2) 
READ(2,.,ENB=M0)  TIB2 
*EA0(2,-)  (TBC(I,2),I=1,NBC) 

CONTINUE 

CALL  INTBC(NX,NBC,DX,TINE,TV) 


SUBROUTINE  CIIKS==========— 

SUBROUTINE  CHKS(NX , NY, I TFS ,S 
DIMENSION  ITFS(NX,2) ,S(NX,2) 


,SO,TV,TVO,TNK,YDIM) 

,S0(NX,2),TV0(NX),tv(- 


>.S(I,2).GE.YDIH)  THEN 


IF(ITFS(I,2). 
ITFS(I,2)  = 0 

end’if  = 0 

teiF(S(I?l)SGT. 
ITFS(I,1)  = 0 

S(  iTiJ’^oT  0 

8(1.2))=  0. 


>.ABS(S(I,2)-S(I,1)).LE. 


IF(Tw(I).GT 
8(1.1) 
SO(I.l) 
8(1.2) 
S0(I,2) 
ITFS(I.l)  = 
ITFS(I,2)  = 


ITFS( I ' 1 ) 


AND. S(  !,!).< 


'.S(I,2).AND.S(I,S 


*S(I.1).AND.S(I,1 


SUBROUTINE  FINDS 


SUBROUTINE  FINDlS(NY,t 
DIMENSION  YI( I) ,TI ( I ) . 
COMMON  /AUL/XDIM.YDIM, 


»SDT  = ABS(S-SO) 

SSO  = S 

IF(DSDT.LE.l.E-08)  THEN 
ODD  = DMAXI(DSDT.  0.0001D0) 
AA(J)°=Jo!5NV 


CC(J) 


CONTINUE 

DD(1)  = 3..(TI(2)-TI(1))/DV 
D0(200  =.23'*<T1(NV-')-t>(Nv))/dv 

?ONnNi^^AA‘J)*(T,(J+1)‘TI(J))tM<J)*(TI(J)-T1(J-,)»/DY 
CALL  TRID(NV,AA,BB,CC,DD) 


if(kt!(;t.ydIm)*rt=ydim  ° 

CONTINUE 

CALL  SPLINE(NY,yi,TI,DD,YL,TI) 
CALL  GREFT(YL,SSO,YPRE,GINT) 


CALL  SPLINE(NY,YI ,TI ,DD,RT,T1 ) 


CALL  SPLINE(NY.YI,TI,DD,RT,T1) 
CALL  GREFT(YLOC,SSO.YPRE,GINT) 


IF(ABS(F).LT™PST)  GOTO  40 
RTO  = (£-FL)/(RT-YL) 
Pif(f.l/o.)  THEN 


IF(((RT-YH).DF-F).((RT-YL).DF-F).GB.O. 

' ■®®^*BS<2-*|f)-6T-*BS(D0.DF))  THEN 

D = O.S.(YH-YL) 

RT  = 0.5*(YH+YL) 


CALL  SPLINE(NY,YI,TI,DD,RT,T1) 
CALL  GREFT(YLOC,SSO,YPRE,GINT) 


20  CONTINUE 
40  CONTINUE 


RETURN 


C=====  SUBROUTINE  TRIO 


SUBROUTINE  TRID(M2,B,D,AXI) 
DIMENSION  A(1),B(I),D(I),CI(1) 


C1(I)=  Cl(I)-RR.Cl(I-l) 

CONTINUE 

R*  = l./D(M2) 


SUBROUTINE  CREFT(YLOC,VPO,VP1,GINT) 
COMMON  /ALL/XDIM,VDIM,DX,DV,DT,CK  CP  RHO 
EPS  = I.E-04 

IP(ABS(yL0C).LT-l.E-O8)  THEN 


1.1415927 


===  SUBROUTINE  FINDT2  =: 


SUBROUTI NE  FINI)T2(  NV . TDK . IIFGCP . 
DIMENSION  VI(NV).T1(NV) 

COMMON  /AU/XDIM.VDIM . DX, DY, DT 
REWIND  (40) 

IF(S.OE.YDIM.OR.S.LE.O.)  RETURN 


=====  SUBROUTINE  FINDT1  ============ 

SUBROUTINE  FINDT1  ( NX , NY,  ncy , HCOv'tINF 

DIMENSION  TO(NX , NY) , T(NX |nY) !tM(NX ) 
DIMENSION  A(I ) ,B( I ) ,C(1 ) ,0(1 ) 

COMMON  /ALL/XDIM , YDI M , DX , DY . DT , CK , CP , 


I»1,JJ)-S..T(I,JJ)+T(I-1,JJ)) 


TO(I,JJ)  = D(J) 


200  CONTINUE 
300  CONTINUE 


SUBROUTINE  OUTPUT(NX,NY,TSTOP,Y 
DIMENSION  Y(NX,NY),T(NX,NY) 
COMMON  /ALL/XDIM,YDIM,DX,DY,DT 
COMMON  /IO/IREAD, IOUT, I REST,  I DU. 


VRITE(IOUT.IOOO)  X,Y(I,J 
CONTINUE 
1 F0RMAT(5E16.0) 


SUBROUTINE  FIND2S(NV ,S01 ,S02 ,S1 ,S2 , IIFGCP ,TMK , 

DIMENSION  YI(1),TI(1)U*(1),BB(I),CC(1),DD(1) 
COMMON  /ALL/XDIM,YDIM,DX,DY,DT,CK,CP,RBO,ALPBA,HFG 


IF(ABS(DSDT1 ) .L 


IF(ABS(DSDT2).L 

IF(S02.LE.l.F.-8 


AA(J)  = 0.5 
BB(J)  = 2.0 
CC(J)  = 0.5 
CONTINUE 

DD(1)  = 3-«(TI(2)-TI(I)  )/DY 
DD(NY)  = 3.»(TI(NY)-TI(NY-1))/DY 


DD(.I)  =3.*(AA(J)*{' 
CONTINUE 

CALL  TRID(NY,AA,BB 


1)-TI(J))+CC(J).(TI(J)-TI<J-I)))/DY 


gint2) 


IF(*BS(YLl-RTl).LT.l.E-8)  RT1  = YL1  + 1 E-OS 
if(ABS(YL2-RT2).lt.l.E-8)  RT2  = YL2  * KE-05 


RETURN 


COMMON  /ALI./XDIM.Y: 


s™n  = TERM  11  -TERM21-EPT«(TERM3I-TERM4 1 ) 

SUHgl  = SUMgl  * SUMIl.SIN(BXL) 

SUMI2  = TERM1 2-TERM22-EPT • ( TERM32- TERM42 ) 
SWIg2  = SUMg2  * SUMI2«SIN(BXL) 


CONTINUE 

1F(ERR.CT.1 

VRITE(.,.) 


SUBROUTINE  FINDTT(NY,TI 


TI(J)  = TI(J)  ♦ 

CONTINUE 

FORMAT('IFH.-I) 


PROGRAM 


n "•  -*  1G"  t“d  llVZj'  T"  *962,  "ThC  ” 


'(eH-O^UU^y^Nlw^ork^pP-  m'-Mo"  **  ' 


sfc  rr»v“- 
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